/*
  Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes related to emission of rays.

// $Id$

#ifndef __emission__
#define __emission__

#include <algorithm> 
#include <vector> 
#include "mcrx-types.h"
#include "ray.h"
#include "sphparticle.h"
#include "montecarlo.h"
#include "random.h"
#include "constants.h"
#include "boost/shared_ptr.hpp"
#include "vecops.h"
#include "angular.h"
#include "grid.h"
#include <memory>

namespace mcrx {
  template<typename, typename> class emission;
  template<typename, typename> class isotropic_emission;
  template<typename, template<typename> class, typename> 
  class emission_collection;
  template <typename, typename> class grid_cell_emission;
  template<typename, typename> class pointsource_emission;
  template<typename, typename> class biased_pointsource_emission;
  template<typename, typename> class laser_emission;
  template<typename, typename> class floodlight_emission;
  template<typename, typename> class divergent_floodlight_emission;
  template<typename, typename> class external_illumination;
  template<typename, typename> class particle_emission;
  template<typename, typename> class particles_emission;

  /** Traits class that defines the class for an emission distribution
      for different grids. */
  template<template <typename> class, typename, typename> class emitter_types {};

  /** Specialization of the emitter_types class for adaptive_grid. */
  template<typename> class adaptive_grid;
  template<typename chromatic_policy, 
	   typename rng_policy> class emitter_types<mcrx::adaptive_grid, 
						    chromatic_policy, 
						    rng_policy> {
  public:   
    typedef grid_cell_emission<chromatic_policy, rng_policy > T_emitter;
  };

  /** Specialization of the emitter_types class for arepo_grid. */
  template<typename> class arepo_grid;
  template<typename, typename> class arepo_cell_emission;
  template<typename chromatic_policy, 
	   typename rng_policy> class emitter_types<mcrx::arepo_grid, 
						    chromatic_policy, 
						    rng_policy> {
  public:   
    typedef arepo_cell_emission<chromatic_policy, rng_policy > T_emitter;
  };
  
}


/** Abstract base class for emission objects, which emit rays
    according to some design probability distribution.  Defines the
    interface.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::emission {
protected:

public:
  typedef typename chromatic_policy::T_lambda T_lambda;
  typedef typename chromatic_policy::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;

  virtual ~emission () {};
  /// Exception class, thrown if there is no emission in the grid.
  class no_emission {}; 
  
  /** The emission function, returns a ray and the angular
      distribution of the emitted rays.  The function takes the policy for
      generating random numbers as an argument, because the random number
      generator is thread-local.  */
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy&)=0; 

  virtual T_lambda zero_lambda() const=0;
  virtual std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const=0;
};


/** Abstract base class for all emission classes with isotropic
    angular distribution.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::isotropic_emission : 
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;

protected:
  isotropic_angular_distribution<T_lambda> ad_;

public:
  isotropic_emission(const T_lambda& l): ad_(l) {};
};


/** The emission_collection class manages emission from a collection
    of emitters. The emitters themselves are not contained in this
    class, it only keeps pointers to them. For emission, it picks the
    individual emitters according to their statistical weight and then
    asks that emitter to emit the ray. */
template <typename chromatic_policy,  
	  template<typename> class sampling_policy, 
	  typename rng_policy_type>
class mcrx::emission_collection :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename mcrx::emission<chromatic_policy, rng_policy_type> T_emission;
  typedef typename T_emission::T_lambda T_lambda;
  typedef typename T_emission::T_biaser T_biaser;
  typedef rng_policy_type T_rng_policy;  
  typedef ray<T_lambda> T_ray;
  typedef probability_distribution<T_emission, 
				   sampling_policy, 
				   T_rng_policy> T_distribution;

protected:
	    
  T_distribution distr_;
  /// Default constructor is only for derived classes.
  emission_collection () {};
  
public:
  emission_collection (const std::vector<emission<chromatic_policy, 
		       T_rng_policy>*>& emitters,
		       const array_1& weights) :
    emission<chromatic_policy, rng_policy_type>()
  { distr_.setup(emitters, weights); };

  /** The emit function simply defers to the individual
      emitters after drawing one. */
  virtual ray_distribution_pair<T_lambda> emit (T_biaser b,
						T_rng_policy& p) {
    return distr_.sample(p)->emit(b, p);};
  T_float total_emission_weight () const {
    return distr_.total_weight();};

  /** Returns a reference to the distribution object. */
  const T_distribution& distribution() const {return distr_;};
  
  /** This function asks the first emitter for its zero_lambda, it's
      dumb but easy. */
  T_lambda zero_lambda() const {
    return distr_.pointers().front()->zero_lambda();};
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new emission_collection<chromatic_policy, sampling_policy, rng_policy_type>(*this));};
};


/** Emission class that handles isotropic emission from a finite sized
    particle. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::particle_emission : 
  public mcrx::isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef particle<T_lambda> T_particle;
  typedef rng_policy_type T_rng_policy;

private:
  /** The particle keeps the data about what is being emitted, like
      the position, velocity and emission of the emitter. */
  T_particle p; 

public:
  particle_emission (const vec3d& position, const vec3d& velocity, 
		     T_float size, T_lambda em):
    isotropic_emission<chromatic_policy, rng_policy_type>(em), 
    p (0, position, velocity, size, em) {};
  
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 

  T_lambda zero_lambda() const {return T_lambda(0.*p.data());}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new particle_emission<chromatic_policy, rng_policy_type>(*this));};
};



/** Emission class representing emission from an isotropic point source. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::pointsource_emission : 
  public mcrx::isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  /// Position of the point source.
  vec3d position;

  T_lambda emission_; ///< The "thing" being emitted.
public:
  pointsource_emission (const vec3d& pos, const T_lambda em):
    isotropic_emission<chromatic_policy, rng_policy_type>(em),
    position (pos), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 

  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new pointsource_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** An isotropic point source which is biased towards the xy
    plane. This does not inherit from isotropic_emission because the
    probability distribution is different, even though it IS
    isotropic. Unfortunately we don't have a nice way of handling
    sources with biased emission, so this is a one-off case. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::biased_pointsource_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

private:
  vec3d position; ///< The position of the source.
  T_lambda emission_; ///< The "thing" being emitted.

public:
  biased_pointsource_emission (const vec3d& pos, 
			       const T_lambda  em ):
    position (pos), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };
  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new biased_pointsource_emission<chromatic_policy, rng_policy_type>(*this));};
};  


/** Emission class representing collimated emission from a point.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::laser_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  /// The angular distribution object.
  collimated_angular_distribution<T_lambda> d;
  
  T_lambda emission_; ///< The "thing" being emitted.

public:
  laser_emission (const vec3d& pos, const vec3d& dir, 
		  const T_lambda  em ):
    position (pos), d (dir, em), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new laser_emission<chromatic_policy, rng_policy_type>(*this));};
};




/** Emission class representing a finite-sized beam of collimated
    emission from a point, covering a certain areal cross-section
    (like an old aircraft searchlight).  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::floodlight_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  T_float radius; ///< The radius of the beam.
  collimated_angular_distribution<T_lambda> d;
  
  T_lambda emission_; // The "thing" being emitted.

public:
  floodlight_emission (const vec3d& pos, const vec3d& dir, T_float r, 
		       const T_lambda  em ):
    position (pos), d (dir, em), radius (r), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new floodlight_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** Emission class representing an almost-collimated beam with
    specified (initial) cross section and divergence. */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::divergent_floodlight_emission :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  vec3d position; ///< The position of the source.
  cone_angular_distribution<T_lambda> d_;
  T_float radius_; ///< The radius of the beam.
  T_float theta_; ///< The opening angle of the radiation.
  T_lambda emission_; // The "thing" being emitted.

public:
  divergent_floodlight_emission (const vec3d& pos, const vec3d& dir, T_float r, 
				 T_float theta, const T_lambda  em ):
    position (pos), d_ (dir, theta, em), radius_ (r), 
    theta_(theta), emission_ (em) {
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new divergent_floodlight_emission<chromatic_policy, rng_policy_type>(*this));};
};


/** Emission class representing an isotropic external illumination
    field, originating on a sphere and directed inwards.  */
template <typename chromatic_policy, typename rng_policy_type>
class mcrx::external_illumination :
  public mcrx::emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;
  
private:
  T_float radius_; ///< The radius of the sphere.
  /** The opening angle of the radiation around the inward direction
      (calculated in the constructor). */
  T_float theta_;
  T_lambda emission_; // The "thing" being emitted.

public:
  /** Creates an external illumination field around the origin. */
  external_illumination (/// Radius where emission originates
			 T_float r, 
			 /** The radius around the origin that should
			     be covered by the emission. */
			 T_float x, 
			 /// The thing being emitted.
			 const T_lambda  em ):
    radius_ (r), theta_(asin(x/r)), emission_ (em) {
    assert(r>=x);
    ASSERT_ALL(em==em);
    ASSERT_ALL(em<1e300);
  };

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new external_illumination<chromatic_policy, rng_policy_type>(*this));};
};


/** Emission class representing isotropic emission from a rectangular,
    axis-aligned grid cell. IMPORTANT: because this object needs a
    pointer to its cell, we can not simply unify cells and put the
    result in the supercell, as the factory does. After that point,
    the cell_ pointer MUST be restored to point to the current
    cell. Neither the factory nor this object can correct this. It's a
    flaw in the inherent assumption that cell data can always be
    independent of its cell, but this is not true when the position is
    needed.  */
template <typename chromatic_policy, 
	  typename rng_policy_type>
class mcrx::grid_cell_emission :
  public isotropic_emission<chromatic_policy, rng_policy_type> {
public:
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_lambda T_lambda;
  typedef typename isotropic_emission<chromatic_policy, rng_policy_type>::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef rng_policy_type T_rng_policy;

  /** To avoid circularly dependent template parameters, it is
  hardcoded that the grid cells contain this class.  This precludes
  making grids which contain objects which in turn contains
  grid_cell_emission objects, but was chosen for simplicity and the
  lack of the current need for anything more advanced. */
  typedef grid_cell<grid_cell_emission> T_cell;

private:
  const T_cell* cell_; ///< Pointer to the cell doing the emission.

  T_lambda emission_; ///< The "thing" being emitted.

public:
  /** Constructor makes an emitter with the specified emission in the
      specified cell. Note that emission_ is a reference copy, because
      the ir_grid relies on the individual emitters referencing the
      main L_lambda array. If this is unacceptable, you need to
      explicitly make an independent_copy before passing it to the
      constructor. */
  grid_cell_emission (const T_lambda& em, const T_cell& c):
    isotropic_emission<chromatic_policy, rng_policy_type>(em),
    cell_(&c), emission_(em) {
    // the value here may be undefined, so we can't check it yet.
  };

  // copy constructor is also needed by factory
  grid_cell_emission(const grid_cell_emission& rhs) :
    isotropic_emission<chromatic_policy, rng_policy_type>(rhs), 
    cell_(rhs.cell_), 
    emission_(independent_copy(rhs.emission_)) {};

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy& p); 
  T_lambda zero_lambda() const {return T_lambda(0.*emission_);}; 
  std::auto_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return std::auto_ptr<emission<chromatic_policy, rng_policy_type> >(new grid_cell_emission<chromatic_policy, rng_policy_type>(*this));};

  /** Returns the cell pointer. This is needed because if cells are
      unified, we need to be able to restore the cell pointers (very
      hoaky) */
  const T_cell* cell() const { return cell_; };
  /** Sets the cell pointer. This is needed because if cells are
      unified, we need to be able to restore the cell pointers (very
      hoaky) */
  void set_cell(const T_cell* c) { cell_= c; };

  // copy constructor and assignment operator are implicitly generated

  const T_lambda& get_emission() const { return emission_; };
  T_lambda& get_emission() { return emission_; };

  /** Post_setup is a no-op for this class. (It is necessary because
      of the arepo_cell_emission class.) */
  static void post_setup(const adaptive_grid<grid_cell_emission>& g) {};

  /// \name grid_factory stuff. 
  /// These are needed for the grid_factory to work.
  ///@{
  grid_cell_emission& operator+= (const grid_cell_emission& rhs) {
    emission_+= rhs.emission_; return *this;};
  grid_cell_emission& operator*= (const grid_cell_emission& rhs) {
    emission_*= rhs.emission_;//content is irrelevant for this operation
    return *this;};
  grid_cell_emission& operator/= (T_float rhs) {
    emission_/= rhs;
    return *this;};
  grid_cell_emission operator*(const grid_cell_emission& rhs) const {
    return grid_cell_emission (*this)*= rhs;};
  grid_cell_emission operator/(T_float rhs) const {
    return grid_cell_emission (*this)/= rhs;};
  static grid_cell_emission unification (const grid_cell_emission& sum, int n) {
    // emission is extensive
    return sum/n;};
  ///@}
};

  

// ***
// *** function definitions ***
// ***


/** Emits a ray from the particle.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::particle_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::particle_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  ASSERT_ALL(p.data()==p.data());
  ASSERT_ALL(p.data()<1e300);

  // get five random numbers
  blitz::TinyVector<T_float, 5> r;
  rng.rnd(r);
  
  // draw a radius based on the SPH smoothing kernel
  const T_float radius = p.draw_radius(r[0]);

  // draw the position in spherical coordinates
  const T_float phi = r [1]*2*constants::pi; 
  const T_float cos_theta = 2*r [2] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [3]*2*constants::pi; 
  const T_float dcos_theta = 2*r [4] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);
  const vec3d raydir (dsin_theta*cos (dphi), 
		      dsin_theta*sin (dphi), 
		      dcos_theta);
  //and create the ray
  boost::shared_ptr<T_ray> rp
    (new T_ray (p.position() + radius*
		vec3d (sin_theta*cos (phi), sin_theta*sin (phi), cos_theta),
		raydir,
		initialize(p.data (), 1.))); // set ray intensity to 0.
  // we do *not* set the doppler factor here because it needs to be
  // emerged to the cameras, etc. instead we return the velocity
  // vector.
  return ray_distribution_pair<T_lambda> (rp, p.data(), this->ad_, p.velocity());
}


/** Isotropically emits a ray from the point source in a random
    direction. The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::pointsource_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::pointsource_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  ASSERT_ALL(emission_==emission_);
  ASSERT_ALL(emission_<1e300);

  // get two random numbers
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  //draw the emission direction (from an isotropic distribution)
  const T_float phi = r [0]*2*constants::pi;  
  const T_float cos_theta = 2*r [1] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);

  // and create the ray
  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (position,
		   vec3d (sin_theta*cos (phi),
			  sin_theta*sin (phi),
			  cos_theta),
		   initialize(emission_, 1.))),
      emission_,
      this->ad_);
}


template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::biased_pointsource_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::biased_pointsource_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  // get 2 random numbers
  blitz::TinyVector<T_float, 2> r;
  rng(r);

  //draw the emission direction (NOT from an isotropic distribution)
  const T_float phi = r [0]*2*constants::pi;  
  // this is a P(cos theta) = pi/4*cos([cos theta]*pi/2) distro
  const T_float cos_theta = 2/constants::pi*asin(2*r [1] - 1);
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);
  //const T_float bias = 0.5*pi*cos(cos_theta*0.5*pi);
  const T_float bias = 2/(constants::pi*cos(cos_theta*0.5*constants::pi));

  // create the ray
  boost::shared_ptr<T_ray> ray (new  T_ray (position,
					    vec3d (sin_theta*cos (phi),
						   sin_theta*sin (phi),
						   cos_theta),
					    initialize(emission_, 1.)));
  // the normalization is bias*emission
  return ray_distribution_pair<T_lambda> (ray, bias*emission_, 
					  this->isotropic_ad);
}


/** Emits a ray from the specified point in the specified
    direction. This is trivial because position and direction are
    already specified, all that is needed is to create the ray.  The
    biaser argument is not used.  */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::laser_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::laser_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (position, d.direction (), initialize(emission_, 1.))), 
       emission_, d);
}


/** Emits a ray from in a beam of specified radius in the specified
    direction.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::floodlight_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::floodlight_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  // get 2 random numbers
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  // We need to draw a position on the emitting disk. The radius is
  // drawn from a square-root distribution.
  const T_float rr = sqrt(r [0])*radius;  
  const T_float phi = r [1]*2*constants::pi;
  // Make a position in the xy plane
  const vec3d xypos(rr*cos(phi), rr*sin(phi), 0.);
  // and rotate+translate so the disk is normal to the desired
  // direction and in the desired position
  const vec3d pos = rotate(xypos, d.direction()) + position;

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos, d.direction (), initialize(emission_, 1.) )), 
       emission_, d);
}


/** Emits a ray from in a beam of specified radius in the specified
    direction.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::divergent_floodlight_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::divergent_floodlight_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  // get 4 random numbers
  blitz::TinyVector<T_float, 4> r;
  rng.rnd(r);

  // We need to draw a position on the emitting disk. The radius is
  // drawn from a square-root distribution.
  const T_float rr = sqrt(r [0])*radius_;  
  const T_float phi = r [1]*2*constants::pi;
  // Make a position in the xy plane
  const vec3d xypos(rr*cos(phi), rr*sin(phi), 0.);
  // and rotate+translate so the disk is normal to the desired
  // direction and in the desired position
  const vec3d pos = rotate(xypos, d_.direction()) + position;

  //draw the emission direction from an isotropic distribution of cos
  //thetas up to theta_, centered around 1 (ie, pointing forward), and
  //rotate into the system of the emission point

  const T_float dphi = r [2]*2*constants::pi;  
  const T_float dcos_theta = 1-(1-cos(theta_))*r [3];
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // Make a vector of the emission direction and rotate it into the
  // correct direction
  vec3d dir(dsin_theta*cos (dphi),
	    dsin_theta*sin (dphi),
	    dcos_theta);
  dir=rotate(dir, d_.direction());

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos, dir, initialize(emission_, 1.) )), 
       emission_, d_);
}


/** Emits a ray in an isotropic direction pointing inward from the
    surface of a sphere.  The biaser argument is not used. */
template <typename chromatic_policy, typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::external_illumination<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::external_illumination<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  // get 4 random numbers
  blitz::TinyVector<T_float, 4> r;
  rng.rnd(r);

  //draw the emission point from an isotropic direction
  const T_float phi = r [0]*2*constants::pi;  
  const T_float cos_theta = 2*r [1] - 1;
  const T_float sin_theta = sqrt (1- cos_theta*cos_theta);

  // Make a vector pointing to the emission point
  const vec3d pos(sin_theta*cos (phi),
		  sin_theta*sin (phi),
		  cos_theta);

  // draw the emission direction from an isotropic distribution of cos
  // thetas up to theta_, centered around -1 (ie, pointing inward),

  const T_float dphi = r [2]*2*constants::pi;  
  const T_float dcos_theta = (1-cos(theta_))*r [3] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // Make a vector of the emission direction and rotate it into the
  // correct coordinate system
  vec3d dir(dsin_theta*cos (dphi),
	    dsin_theta*sin (dphi),
	    dcos_theta);
  dir=rotate(dir, pos);

  // generate a cone angular distribution
  const cone_angular_distribution<T_lambda> d(dir, theta_, emission_);

  return ray_distribution_pair<T_lambda> 
    ( boost::shared_ptr<T_ray> 
      (new  T_ray (pos*radius_, dir, initialize(emission_, 1.) )), 
       emission_, d);
}


/** Emits a ray from the grid cell. The biaser argument is not used. */
template <typename chromatic_policy,
	  typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::grid_cell_emission<chromatic_policy, rng_policy_type>::T_lambda> 
mcrx::grid_cell_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng)
{
  ASSERT_ALL(emission_==emission_);
  ASSERT_ALL(emission_<1e300);

  // get five random numbers
  blitz::TinyVector<T_float, 5> r;
  rng.rnd(r);
  
  // randomly draw a position
  const vec3d p (r [0], r [1], r [2]);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [3]*2*constants::pi; 
  const T_float dcos_theta = 2*r [4] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  //and create the ray
  return ray_distribution_pair<T_lambda> 
    (boost::shared_ptr<T_ray> 
     (new T_ray (cell_->getmin() + p*cell_->getsize(),
		 vec3d (dsin_theta*cos (dphi), dsin_theta*sin (dphi), 
			dcos_theta),
		 initialize(emission_, 1.))), 
     emission_,
     this->ad_);
}

#endif
